Arteriovenous malformation Map2k1 mutation affects vasculogenesis

Somatic activating MAP2K1 mutations in endothelial cells (ECs) cause extracranial arteriovenous malformation (AVM). We previously reported the generation of a mouse line allowing inducible expression of constitutively active MAP2K1 (p.K57N) from the Rosa locus (R26GT-Map2k1-GFP/+) and showed, using Tg-Cdh5CreER, that EC expression of mutant MAP2K1 is sufficient for the development of vascular malformations in the brain, ear, and intestines. To gain further insight into the mechanism by which mutant MAP2K1 drives AVM development, we induced MAP2K1 (p.K57N) expression in ECs of postnatal-day-1 pups (P1) and investigated the changes in gene expression in P9 brain ECs by RNA-seq. We found that over-expression of MAP2K1 altered the transcript abundance of > 1600 genes. Several genes had > 20-fold changes between MAP2K1 expressing and wild-type ECs; the highest were Col15a1 (39-fold) and Itgb3 (24-fold). Increased expression of COL15A1 in R26GT-Map2k1-GFP/+; Tg-Cdh5CreER+/− brain ECs was validated by immunostaining. Ontology showed that differentially expressed genes were involved in processes important for vasculogenesis (e.g., cell migration, adhesion, extracellular matrix organization, tube formation, angiogenesis). Understanding how these genes and pathways contribute to AVM formation will help identify targets for therapeutic intervention.

Principal component analysis grouped the R26 GT-Map2k1-GFP/+ ; Tg-Cdh5CreER +/− RNA seq datasets separately from the wild-type datasets (Fig. 2). Because EC cell recovery, RNA extraction, library preparation, and sequencing were performed on all samples at the same time, the difference between mutant and wild-type EC datasets likely represents biologic differences rather than artifactual batch effects. Differential gene expression analysis identified 790 genes whose transcript abundance increased, and 887 genes whose transcript abundance decreased (Fig. 3A and Supplementary Tables 1, 2). Map2k1 was identified among the increased transcripts. To estimate the relative expression of the recombined R26 GT-Map2k1-GFP/+ allele compared to the two endogenous alleles, we counted mutant versus wild type reads and normalized for the % of recombined cells. This analysis determined that mRNA expression from a recombined R26 GT-Map2k1-GFP/+ allele was on average 69-fold higher than endogenous Map2k1 expression (Fig. 3B).

Discussion
Previous AVM RNA sequencing studies have involved either human brain AVM specimens which do not contain MAP2K1 mutations or in vitro engineered KRAS (p.G12D) overexpressing ECs [14][15][16][17] . In these systems genes involved in inflammation, cell migration, extracellular matrix remodeling, and angiogenesis are differentially expressed. Our study illustrates the transcriptional profiling of ECs from MAP2K1 (p.K57N) induced brain vascular malformations using a R26 GT-Map2k1-GFP/+ ; Tg-Cdh5CreER +/− mouse line. Our RNA sequencing results are consistent with single cell RNA data from human brain AVMs which also exhibited increased COL15A1, PGF, PLVAP, STC1, ANGPT2 and decreased MFSD2A, SLC16A1, SLC38A5 16,17 . Advantages of our model are that it is in-vivo, has a known mutation, and produces an obvious phenotype. Disadvantages of our mouse line include the possibility that gene expression changes might be different if mutant MAP2K1 was expressed at endogenous levels or if changes were compared to mice expressing wild type MAP2K1 protein from the Rosa locus. The transcriptome changes in MAP2K1 expressing ECs provides insights into the mechanisms that contribute to AVM formation. Overexpression of MAP2K1 resulted in more than 1600 genes whose expression was significantly different between mutant and non-mutant ECs. None of the genes alone appear to be necessary and sufficient to produce an AVM. However, many of the overexpressed transcripts in our dataset promote vasculogenesis. For example, upregulated Esm1 and Vav3 expression may cause dysregulated EC migration and filopodia extension. ECs in Esm1 null mice exhibit delayed vascular outgrowth and reduced filopodia extension 18 . Microvascular ECs overexpressing VAV3 demonstrate increased migration and assembly while embryonic fibroblasts deficient in VAV3 have reduced filopodia extension 19 . COL15A1 forms a proangiogenic scaffold for ECs. ECs cultured with fibroblasts overexpressing COL15A1 have increased migration and tube formation 20 . ITGB3 is a heterodimeric adhesion receptor that interacts with the extracellular matrix and is highly expressed in vessels undergoing pathological angiogenesis 21 . ROBO1, CDH13, PGF, ANGPT2, and FGFR1 increase EC proliferation, migration, and tube formation [22][23][24][25] . Dysregulated expression of multiple proangiogenic factors in MAP2K1 expressing ECs may result in their inability to form normal capillary beds which could lead to arteriovenous shunts.
It is not feasible to obtain mice with a global or conditional knockout allele for every differentially expressed gene to test that gene's role in producing vascular malformations when mutant Map2k1 is expressed. However, it is possible to use genes with both large fold-changes in transcript abundance and increased protein expression for cell-based assays of mutant Map2k1 effects. For example, CD109, whose transcript is 13-fold more abundant in mutant ECs, can be used to flow sort ECs in perturb-seq experiments. As an alternative to antibody-mediated flow sorting, a fluorescent reporter could be knocked into the Col15a1 locus, whose transcript is 39-fold more abundant in MAP2K1 mutant ECs.
Our results showing that the Map2k1 mutation in ECs affects vasculogenesis genes and pathways contributes to our understanding of AVMs pathogenesis. The activating MAP2K1 EC mutation may cause AVM enlargement through abnormal blood vessel production that secondarily causes overgrowth of tissues. This would explain the finding that human AVMs enlarge by a cell non-autonomous mechanism 26 . Drugs that inhibit EC migration, adhesion, and angiogenesis may prove beneficial to patients with AVMs. Anti-angiogenic medications might prevent enlargement of AVMs as well as reduce their recurrence following embolization and resection. Combining anti-angiogenic drugs with MAP2K1 inhibitors for AVM may prove synergistic, improving the efficacy of drugs like trametinib and reducing its toxicity.

Histology and immunostaining.
Endothelial cell isolation. ECs were isolated as previously described 6 . Briefly, whole brain was minced with a scalpel, digested with collagenase A (Roche # 11088793001) at 37 °C for 1.5 h, and passed through an 18-gauge needle. Cells were filtered through a 70-μm strainer and incubated for 15 min at 4 °C with anti-CD31 antibody (BD Biosciences #553370) that had been previously conjugated to magnetic dynabeads (ThermoFisher #11035) following the manufacturer's protocol. Non-bound cells were removed by washing three times with 0.1% BSA in PBS for one minute/wash at room temperature. The bound cells were resuspended in PBS and 1/20th of the suspension was kept for DNA extraction.
Droplet digital PCR (ddPCR). The R26 GT-Map2k1-GFP allele contains from 5′ to 3′: a CAG promoter, a LoxP flanked stop cassette (gene trap), the Map2k1 (p.K57N) cDNA, an IRES-GFP cassette and a poly A signal 13 . DNA was extracted from isolated brain ECs using the DNeasy blood and tissue kit (Qiagen #69506). To detect the non-recombined R26 GT-Map2k1-GFP allele, a ddPCR assay was designed using a forward primer (FP) immediately upstream of the 5′ LoxP side of gene trap (5′-cgtgctggttattgtgctgtc-3′) and a reverse primer (RP) inside the gene trap (5′-gcgttggctacccgtgata-3′). The resulting amplicon size was 278 bp. A gene trap specific probe (5hex/caacac ccgtgcgttttattctgt/3IABkFQ) was contained within the amplicon. To detect the recombined R26 GT-Map2k1-GFP allele a ddPCR assay was designed using the same FP in combination with a RP located within the Map2k1 (p.K57N) cDNA (5′-gatctccaggtggatcagctt-3′ . Reads were normalized using DESeq2's median of ratios, a sample-specific normalization method determined by the median ratio of gene counts relative to geometric mean per gene. Principal component analysis was performed to assess similarities and differences across the samples. Differential gene expression analysis was carried out using DESeq2. We focused our attention on the top 80% of these ranked by average normalized read counts across all eight samples. RNA sequencing data are available in the NCBI gene expression omnibus (GEO) data repository.
Statistical methods. Significantly differentially expressed genes had adjusted p values < 0.05 (after controlling for multiple comparisons using the Benjamani-Hochberg Method, and mean log2-fold change > or < 1). Pathway analysis using the top 200 differentially expressed genes (ranked by adjusted p value) was performed using the Gene Ontology Resource 29 .

Data availability
The datasets generated and/or analyzed during the current study has been submitted to the Gene Expression Omnibus repository, accession number: GSE226484 (https:// www. ncbi. nlm. nih. gov/ geo/ query/ acc. cgi? acc= GSE22 6484).